devtools::install("../../elasticproc2d", quiet=TRUE, quick=TRUE)
library(elasticproc2d)
library(ggplot2)
library(grid)
library(gridExtra)
library(dplyr)
library(cowplot)
# Display plots in smaller size
options(repr.plot.width=15, repr.plot.height=5)
Load some datasets with random rotation and scaling.
source("../scripts/load_datasets.R")
set.seed(18)
spirals <- curves.spiral(n_curves=10, rotate=TRUE, scale=TRUE, center=TRUE)
digits3 <- curves.digit3(rotate=TRUE, scale=TRUE, center=TRUE)
Normalize curves.
Take a look at the data curves.
p1 <- ggplot(bind_rows(spirals, .id="id"), aes(x=X1, y=X2)) +
geom_path(size=0.5, aes(group=id), color="grey") +
coord_fixed()
p2 <- ggplot(bind_rows(digits3, .id="id"), aes(x=X1, y=X2)) +
geom_path(size=0.5, aes(group=id), color="grey") +
coord_fixed()
grid.arrange(p1, p2, nrow=1, top="Spirals and digits3 with random rotation and scaling")
Define functions for plotting.
plot_cov_surface<- function(cov_fit, h=0.01){
# Define covariance surface grid (s,t).
arg.grid = seq(0, 1, by=h)
cov.grid = expand.grid(t = arg.grid, s = arg.grid)
# Evaluate fit on grid.
cov.re = predict(cov_fit$re, newdata = cov.grid)
cov.im = predict(cov_fit$im, newdata = cov.grid)
par(mfrow=c(1,2), mar=c(4,4,4,1), oma=c(0.5,0.5,0.5,0))
# From 'fdapace/src/R/CreateCovPlot.R'
args1 <- list(
xlab='t', ylab='s', zlab = 'C(s,t)',
main = 'Smoothed covariance surface (real part)'
)
args2 = list (x = arg.grid, y = arg.grid, z = matrix(cov.re, nrow=101))
do.call(plot3D::persp3D, c(args2, args1))
# From 'fdapace/src/R/CreateCovPlot.R'
args1 <- list(
xlab='t', ylab='s', zlab = 'C(s,t)',
main = 'Smoothed covariance surface (imaginary part)'
)
args2 = list (x = arg.grid, y = arg.grid, z = matrix(cov.im,nrow=101))
do.call(plot3D::persp3D, c(args2, args1))
}
library(grid)
plot_all <- function(mean, h=0.01){
cov_fit <- mean$fit$cov_fit
# Define covariance surface grid (s,t).
arg.grid = seq(0, 1, by=h)
cov.grid = expand.grid(t = arg.grid, s = arg.grid)
# Evaluate fit on grid.
cov.re = predict(cov_fit$re, newdata = cov.grid)
cov.im = predict(cov_fit$im, newdata = cov.grid)
par(mfrow=c(1,3), mar=c(4,4,4,1), oma=c(0.33,0.33,0.33,0))
# From 'fdapace/src/R/CreateCovPlot.R'
args1 <- list(
xlab='t', ylab='s', zlab = 'C(s,t)',
main = 'Smoothed covariance surface (real part)'
)
args2 = list (x = arg.grid, y = arg.grid, z = matrix(cov.re, nrow=101))
do.call(plot3D::persp3D, c(args2, args1))
# From 'fdapace/src/R/CreateCovPlot.R'
args1 <- list(
xlab='t', ylab='s', zlab = 'C(s,t)',
main = 'Smoothed covariance surface (imaginary part)'
)
args2 = list (x = arg.grid, y = arg.grid, z = matrix(cov.im,nrow=101))
do.call(plot3D::persp3D, c(args2, args1))
plot(mean)
}
order <- list("polygon", "smooth")
penalty <- list(0, 1, 2)
# Turn of warnings
defaultW <- getOption("warn")
options(warn = -1)
spiral.means <- lapply(order, function(type){
lapply(penalty, function(pen){
start_time <- Sys.time()
knots <- seq(0, 1, length = if_else(type == "smooth", 13, 19))
print(paste0("Computing spirals mean: Type = ", type, " Penalty = ", pen, " #Knots = ", length(knots)))
flush.console()
mean <- compute_elastic_proc2d_mean(spirals, knots = knots, type = type, penalty = pen)
print(Sys.time() - start_time)
flush.console()
mean
})
})
[1] "Computing spirals mean: Type = polygon Penalty = 0 #Knots = 19" Time difference of 1.163191 mins [1] "Computing spirals mean: Type = polygon Penalty = 1 #Knots = 19" Time difference of 44.74742 secs [1] "Computing spirals mean: Type = polygon Penalty = 2 #Knots = 19" Time difference of 20.94888 secs [1] "Computing spirals mean: Type = smooth Penalty = 0 #Knots = 13" Time difference of 11.57825 secs [1] "Computing spirals mean: Type = smooth Penalty = 1 #Knots = 13" Time difference of 11.86377 secs [1] "Computing spirals mean: Type = smooth Penalty = 2 #Knots = 13" Time difference of 11.56085 secs
void <- lapply(1:length(order), function(i){
lapply(1:length(penalty), function(j){
mean <- spiral.means[[i]][[j]]
plot_all(mean)
print(paste0("Type = ", mean$type, " Penalty = ", penalty[[j]], " #Knots = ", length(mean$knots)))
})
})
[1] "Type = polygon Penalty = 0 #Knots = 19"
[1] "Type = polygon Penalty = 1 #Knots = 19"
[1] "Type = polygon Penalty = 2 #Knots = 19"
[1] "Type = smooth Penalty = 0 #Knots = 13"
[1] "Type = smooth Penalty = 1 #Knots = 13"
[1] "Type = smooth Penalty = 2 #Knots = 13"
digits3.means <- lapply(order, function(type){
lapply(penalty, function(pen){
start_time <- Sys.time()
knots <- seq(0, 1, length = if_else(type == "smooth", 13, 19))
print(paste0("Computing digits3 mean: Type = ", type, " Penalty = ", pen, " #Knots = ", length(knots)))
flush.console()
mean <- compute_elastic_proc2d_mean(digits3, knots = knots, type = type, penalty = pen)
print(Sys.time() - start_time)
flush.console()
mean
})
})
[1] "Computing digits3 mean: Type = polygon Penalty = 0 #Knots = 19" Time difference of 56.27127 secs [1] "Computing digits3 mean: Type = polygon Penalty = 1 #Knots = 19" Time difference of 20.34927 secs [1] "Computing digits3 mean: Type = polygon Penalty = 2 #Knots = 19" Time difference of 20.33894 secs [1] "Computing digits3 mean: Type = smooth Penalty = 0 #Knots = 13" Time difference of 32.41793 secs [1] "Computing digits3 mean: Type = smooth Penalty = 1 #Knots = 13" Time difference of 32.61988 secs [1] "Computing digits3 mean: Type = smooth Penalty = 2 #Knots = 13" Time difference of 32.37099 secs
void <- lapply(1:length(order), function(i){
lapply(1:length(penalty), function(j){
mean <- digits3.means[[i]][[j]]
plot_all(mean)
print(paste0("Type = ", mean$type, " Penalty = ", penalty[[j]], " #Knots = ", length(mean$knots)))
})
})
[1] "Type = polygon Penalty = 0 #Knots = 19"
[1] "Type = polygon Penalty = 1 #Knots = 19"
[1] "Type = polygon Penalty = 2 #Knots = 19"
[1] "Type = smooth Penalty = 0 #Knots = 13"
[1] "Type = smooth Penalty = 1 #Knots = 13"
[1] "Type = smooth Penalty = 2 #Knots = 13"
mean.sp.nofx.poly <- compute_elastic_proc2d_mean(spirals, knots = seq(0,1,length=11), type = "polygon", penalty = -1)
mean.sp.nofx.smooth <- compute_elastic_proc2d_mean(spirals, knots = seq(0,1,length=11), type = "smooth", penalty = -1)
mean.d3.nofx.poly <- compute_elastic_proc2d_mean(digits3, knots = seq(0,1,length=11), type = "polygon", penalty = -1)
mean.d3.nofx.smooth <- compute_elastic_proc2d_mean(digits3, knots = seq(0,1,length=11), type = "smooth", penalty = -1)
Plots
plot_all(mean.sp.nofx.poly)
plot_all(mean.sp.nofx.smooth)
plot_all(mean.d3.nofx.poly)
plot_all(mean.d3.nofx.smooth)